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ABSTRACT 


An  algorithm  has  been  developed  for  finding  the  eigenvalues  of  a 
symmetric  matrix  A  in  a  given  interval  [a,  b]  and  the  corresponding 
eigenvectors  using  a  modification  of  the  method  of  simultaneous  itera- 
tions with  the  same  favorable  convergence  properties.   The  technique 
is  most  suitable  for  large  sparse  matrices  and  can  be  effectively  im- 
plemented on  a  parallel  computer  such  as  the  ILLIAC  IV. 
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Introduction 

Consider  the  eigenvalue  protlem  Au  =  Xu  where  A  is  a  real  symmetric 

matrix  of  order  n.   If  n  is  large  and  A  is  sparse,  then  an  iterative 

method  for  finding  the  eigenvalues  and  eigenvectors  which  uses  A  only  as 

an  operator  may  be  superior  to  a  method  which  reduces  A  to  a  condensed 

form  (e.g.,  tridiagonal)  using  orthogonal  transformations  resulting  in  a 

far  denser  matrix  exceeding  storage  capacity.   This  superiority  is  quite 

demonstrable  if  only  a  few  of  the  eigenvalues  and  the  corresponding 

eigenvectors  are  required.   Bauer's  method  of  Simultaneous  Iterations 

[l,  2]  for  finding  a  few  of  the  leading  eigenvalues  and  eigenvectors  is 

one  of  such  methods.   In  this  paper  we  modify  his  method  to  find  the 

eigenvalues  of  A  in  some  interval  [a,  b].   Without  loss  of  generality  we 

assume  that  A  is  positive-definite.   The  problem  is  therefore  to  find  the 

q  eigenvalues  X,Ajt,...,X,   ^,  where 
p   p+1        p+q-1 

c„>X  >A^>  . . .  >X   >b>X  >X  ^  >  ...  >X  ,    >a>X  .  >  ...  >X  >c  >0 
2—  1—  2—     —  p-1—   p—  p+1—     —  p+q-l  —  p+q—     —  n—  1 

and  the  corresponding  eigenvectors . 

The  Algorithm 

Similar  to  the  method  of  "Simultaneous  Iterations"  [2],  we  propose 

the  following  iterative  scheme: 

(i)   Choose  an  nxq  matrix  X  such  that 

o 

X^  X  =  I 
o   o    q 

^ii)   \+k  =  V^)\         k=  0,  1,  2,  ...  (1) 

where  B  is  a  polynomial  of  A  to  be  defined  later,  and  T  (B)  is  the 

m 

Chebyshev  polynomial  of  B  of  degree  m. 

(iii)  Using  Gram-Schmidt  orthonormalization  process  we  decompose  X  as 

J 


where  U.U.  =  I  and  R.  is  an  upper  triancxilar  matrix  of  order  q. 
(iv)  Let 

Z^=AU^.  (3) 

Then  form  the  positive-definite  qxq  matrix 

and  solve  for  its  eigenvectors  Q  . 

J 

(v)  Thus, 

Go  back  to  (ii)  and  so  on  until  X  AX  approaches  a  diagonal  matrix  whose 
elements  are  the  eigenvalues  of  A  in  [a,  b];  this  occurs  when  Q  ap- 
proaches  the  identity  matrix. 

The  matrix  B  in  (l)  may  be  obtained  as  follows: 
Let 

A  =  [2A  -  (a+b)l]/(b-a);  (6) 

then  an  eigenvalue  of  A  is  given  by 

X  =  [2A  -  (a+b)]/(b-a)  (?) 

and  those  eigenvalues  of  A  in  [-1,  l]  correspond  to  the  eigenvalues  of 
A  in  [a,  b].   The  interval  [-c,  c]  that  contains  all  the  eigenvalues  of 
A  is  given  by, 

c  =  max{d  ,  d  }  (8) 

where 

^1  "  l^^i  "  (a+b)|/(b-a) 

d^  =  |2c2  -  (a+b)|/(b-a) 
If  some  mapping  y  =  f(x)  can  be  found  such  that  f{\)  is  outside  the  inter- 
val [-1,  l]  for  X  in  [-1,  l]  and  |f(X)|<_l  for  X  outside  [-1,  l],  then  the 
algorithm  described  above  (i)-(v)  will  yield  the  eigenvalues  of  A  in 


(9) 


[a,  b]  and  the  corresponding  eigenvectors. 

Let  us  therefore  map  the  interval  [-1,  l]  on  the  subintervals  [-c,  -l] 
and  [l,  c],  one  such  mapping  may  be  taken  as, 

y  =  ±/ax+3 
or       X  =  (y^-6)/a.  (lO) 

X  =  -1  corresponds  to  y  =  ±1,  thus 

1  =  -a+3; 

and  X  =  +1  corresponds  to  y  =  ±c,  thus 

2 
c  =  a+3. 

Therefore 


and 


(11) 


a  =  |(c2-l)   • 
3  =  |(c2+l) 
Hence  the  matrix  B  is  taken  as 

B  =  [2A^  -  (c^+l)l]/(c^-l)  (12) 


y  =  A(B)  =  [2X^  -  (c^+l)]/(c^-l)  (l3) 


i.e.,     -(^)  <  y  <  -1. 
a 

We  now  introduce  the  following  theorems. 
Theorem  1: 

Let  E  be  the  linear  space  spanning  the  columns  of  X  .   In  case  of 
stable  convergence  (no  reordering  of  the  eigenvalues  if  the  LR-Cholesky  de- 
composition is  applied  to  G  ),  the  angle  4).    between  the  i-th  eigenvector 
u.  and  the  linear  space  E  =  {x|x=T  (B)y,  yeE   },  spanning  the  columns  of 

-L  fj  111  (J  ""-L 

X  5  is  asymptotically  for  j-x»  of  order  0(q,  )  in  which 
J  -^ 

q.  =  max  {  |T   (y  )|/|T   (y  )|}  (ll;) 

1   k?fP    "^   k  I  I  m   1  I 

ieP,  where  P  =  {p,  p+1,  ...,  p+q-l}.   The  proof  is  quite  similar 


to  that  of  Theorem  2  in  [2]  and  hence  will  be  omitted  here  (see  Appendix 

I). 

Theorem  2: 

The  columns  of  the  matrices  X  as  generated  by  (i)-(v)  are  such  that 

l|u.  -  x.^-^^ll  =  O(q.J)  (15) 

where  q.  is  as  given  by  (lU).   The  proof  is  again  similar  to  that  of 
Theorem  (3)  in  [2]  (see  Appendix  l). 

A  proper  order  of  the  Chebyshev  polynomial,  m,  can  be  obtained,  as 

in  [2],  by  stipulating  that  the  parallelization  of  the  columns  of  X, 

k+m 

should  not  go  further  than  that  at  most  one  decimal  digit  is  cancelled 
out  when  these  columns  are  orthonormalized,  i.e., 

l\-i  '^^>l  '^°' 

where 

|y^|  =  max  |y^| .  (l6) 

ieP 

An  approximation  of  y   is  obtained  from  (13)  by  replacing  A  by 

1 

[2A  (G  )  -  (a+b)]/(b-a)  where  A  (G  )  is  the  maximum  eigenvalue  of  G  . 
■■^  J  ^  J  J 

Thus, 

cosh  [(m-l)  cosh    IPpI]  <  10 >  and 

T  ^   cosh"  10   „     3 /,„>, 

m-l  <  — =  — (ITJ 

cosh    |y^|    cosh"   |y^| 
As  soon  as  one  of  the  eigenvalues  of  G  stagnates,  the  corresponding 

J 

eigenvalue  of  A  can  be   found  within  computer  accuracy.      Once  this  happens 
we  can  test   for  acceptance  of  the   corresponding  eigenvector.      The  accep- 
tance test   is   that   of   [3]   except  that  the   discoTinting  rule   is   given  by 

f.    :=  f.    X  /p^tP  (18) 

1  1        T    (a.  ) 

m     1 


where  a.    is  an  approximation  to  y,  obtained  as  explained  above,  and  h  is 
the  number  of  eigenvalues  already  accepted. 

Some  Numerical  Results  . 

The  algorithm  described  above  has  been  implemented  in  Fortran  on 
UCLA's  IBM  360/91  computer.   The  flow  chart  is  slightly  different  from  that 
of  Rutishauser  [3], 


zl:=z2:=ks:=0 
m:=2 


Perform  m  steps  with 

Chebyshev  iterations 

ks: =ks+m 


z2:=z2-m 


Yes 


zl:=zl+l 

z2:=Uxzl 

m:  =2xni 


TRED2,  IMTQL2  [h] 
ks:=ks+l 


i 


Acceptance  tests,  reduce  m  if  needed 


To  demonstrate  the  numerical  behavior  of  the  algorithm  we  tested  three 
symmetric  matrices  A  ,  A  ,  and  A   (shown  in  Appendix  II).   In  all  three 
examples  we  obtained  the  required  eigenvalues  correct  to  lU  decimal  places 
and  the  eigenvectors  correct  to  at  least  7  decimal  places. 
Example  1 

The  Gerschgorin  disks  of  the  Sh'xSk   matrix  A  show  that  we  have  16 
eigenvalues  in  the  interval  (l.i+,  3.6),  8  eigenvalues  in  the  interval 
(5.^,  ^'6)i   and  i+0  eigenvalues  in  the  interval  (8.^,  16.O).  We  seek  to 


obtain  the  8  intermediate  eigenvalues  and  eigenvectors.   Two  different 
intervals  [a,  b]  have  been  tested  [U.O,  8.0],  and  [5-3,  6.7]  with 

^o  "  ^^kV    %2'  •••'  %8^- 

TABLE  I 
[a,  b]     M     zl     ks 

[i+.O,  8.0]   l6     T     lli+ 
[5.3,  6.7]   32    15    U88 

M:   maximum  degree  of  Chebyshev 
polynomials 

zl:  niomber  of  QL  steps 

ks:  nximber  of  iteration  steps 

Table  I  shows  that,  for  the  same  acceptance  test,  the  number  of  itera- 
tion steps  required  (ks )  is  smaller  when  a  and  b  are  as  far  as  possible 

from  the  eigenvalues  to  be  evaluated  (X  ,  A  ,_,  ....  A    .),  which  is 

p   p+1        P+q-1 

clear  from  relation  (l^)  and  the  preceeding  interval  transformation. 

If  we  choose  X  =  [e,^,  e,  ,  ...,  e,  ] ,  however,  and  [a,  b]E[5.3,  6.7] 
the  process  converges  to  the  required  8  eigenvalues  and  eigenvectors  in 
only  265  iteration  steps  (ks  =  265),  zl  =  10,  but  M  =  58. 
Example  2 

The  6'x6   positive-definite  matrix  A  has  the  simple  roots  1,  5,  25, 
and  a  triple  root  at  15.   If  we  seek  the  triple  root  15  and  take  the  in- 
terval [a,  b]  as  [7,  2k]   we  require  62  iteration  steps,  zl  =  5,  and  the 
maximum  degree  M  of  the  Chebyshev  polynomials  is  6.   Here  we  only  obtain  an 
invariant  subspace  that  corresponds  to  X(A  )  =  15 . 

If,  however,  we  seek  the  four  eigenvalues  in  the  interval  [2,  2k] 
we  require  only  20  iteration  steps  with  the  degree  of  Chebyshev  polynomials 
not  exceeding  2,  and  zl  =  3.   Again  we  obtain  an  invariant  subspace  corre- 


spending  to  X(A    )   =  15  and  the  correct   eigenvector  corresponding  to   X(A    )   =  5. 


Example  3 


i+r        kTT 


The  positive-definite  matrix  A     has   the  eigenvalues   X     =  l6  sin   [    /    . -■  \  ] , 
k  =  1,    2,    ...,    n,   where  n  =  6h.      There  are  26  eigenvalues   in    (0,    2),    6   in 
(2,    h)y   and  32  in   (U,   l6).      Here  we  woiold  like  to  evaluate  the  eigenvalues 
in  the   interval   [2,    h]   and  the  corresponding  eigenvectors.      Table  II  shows 
the  effect  of  our  ass\unption  regarding  the  distribution  of  eigenvalues   on 
the  number  of  iteration  steps    (ks )   required  for  convergence  to  the  eigen- 
values  and  eigenvectors   in  the  interval   [2,   i+],    and  degrees  of  the  Chebyshev 
polynomials . 


TABLE 

II 

X 
o 

M 

zl 

ks 

(i)  [633,  .. 

■•'  ^38^ 

52 

9 

290 

(ii)  [e^-|_>  •• 

•-  ^0^ 

52 

8 

237 

(iii)  [e^^,  .. 

••'  "37^ 

52 

21 

1312 

(iv)  [e^^,  .. 

'•'  ^28^ 

52 

10 

3U3 

Experiment    (iii)   indicates   that  using  fewer  columns   in  X     than  the  actual 

o 

number  of  eigenvalues  in  [2,  h]   led  to  a  large  ntimber  of  iteration  steps 
ks  =  1312,  and  zl  =  21  to  converge  to  the  assumed  3  eigenvalues.  While 
experiment  (iv)  shows  that  if  we  use  more  columns  in  X  than  the  actual 
number  of  eigenvalues,  even  if  we  completely  misjudge  the  distribution  of 
the  eigenvalues  in  the  intervals  (O,  2),  (2,  k),    (k^   16),  the  n\amber  of 
iteration  steps  (ks)  and  the  number  of  QL  steps  (zl)  required  for  conver- 
gence to  the  6  eigenvalues  and  eigenvectors  in  [2,  h]   are  only  slightly 
more  than  those  in  experiment  (i)  where  we  used  the  true  distribution  of 


the  eigenvalues.   Since  the  nvimber  of  iterations  required  appears  to  be 
largely  independent  of  the  choice  of  coliomns  of  the  initial  matrix  X  ,  a 
random  number  generator  could  be  used  to  generate  these. 

We  notice  that  the  above  algorithm  is  quite  suitable  for  parallel 
computations  since  the  major  operation  here  is  multiplying  a  matrix  by  a 
vector,  which  can  be  handled  rather  efficiently  on  a  parallel  computer 
such  as  the  ILLIAC  IV. 
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Appendix  I 


Proof  of  Theorem  1 


The  iterations  (i)-(v)  are  orthogonally  invariant,  i.e.,  replacing 

A  by  h''^AH  =  A  (where  h''^H  =  I,  A  =  diag  (A.))  and  X  by  h"'^X  has  the  ef- 
•^  1        o       o 

feet  that  all  X  are  replaced  by  H  X  ,  while  the  G  and  X  are  not 
J  J  J      J 

changed.   Therefore  we  can  ass\ime,  without  loss  of  generality,  that 


A  =  diag  (A^,  A^,  ...,  A^). 
In  case  of  stable  convergence  E  can  be  spanned  by  q  vectors 


^1 


p-1,1 


n,l 


^2' 


"l,. 


p-1,2     p-l,q 


n,2' 


••1 


^p+q,l     .p+q,2    .p+q,q 


n,q 


row  p 


(A.l) 


row  p+q-1 


10 


According  to    (l),    E.    is   spanned  by 

J 


\^^^l)    -11 


m       p-1       p-1,1 


m       p 


T   "^(y   _^    )    X  ^      n 
m    .  p+q        P+q,l 


T  "^(y    )    X     . 
m       n       n,l 


\^[^1^    \2 \^^)    ^,, 


Tj(y         )    x^   ,    _ 
m       p-1        p-l>2 


m       p+1 


J 


T  "^(y     J   X     , 
m       p-1       p-l,q 


m       p+q-1 


T  "^  (y        )   X  ...  T  "^  (y        )   x 

m    .  P+q       P+q, 2  m    .  p+q       p+q,q 


T  "^(y    )    X     ^ T  hv    )   X 

mnn,2  mnn,q 


(A. 2) 


For  example,   the  angle  between  e     and  the   first   column  of   (A. 2)    is   given 

XT 


cos        m        = 
^P 


T  '^hv  ) 


m 


I,  K'K^  ^,.y 


P  =    {p+l,    p+2,    . .. ,   p+q-1} 


hence. 


.    2 


=^"  *p  =  Jp  <V(\)  \,i>^Jp  'Vf\'  %,i'  •       ^  '  'P-  ^' 


or 


where 


'i'     is  of  order  0(q_'^) 
P  T> 

q^  =  max  (  |  T   (y    )  |  /  |  T   (y    )  |  },  k  ?f   P 

^         ^        '    m     k    '    '    m     p    ' 


11 


Therefore, 

<t>.^      is  at  most  of  order  OCq.*^)   | 
Proof  of  Theorem  2 

Taking  the  q  vectors  (A.l)  each  divided  byT"^(vi),  £  =  p,  p+1, 
....  p+q-1,  as  coordinate  vectors  W  =  [w  ,  w  ,^ ,  . . . ,  w  .   .]  in  E,,  the 

p   p+1      P+q-i     J 

-2 
eigendirections  of  the  projected  operator  A   are  the  nxq  matrix  Y  =  WS 


where 


in  which. 


y"^  A"2  Y  =  Dj"^  (A.  3) 


D.  =  diagld    »  d      ,  ...,  d  .   ^    j, 
J        P     P+1         P+q-1 

Y^  Y  =  I  , 

q 


and  S   is   a  qxq  orthogonal  matrix.      Thus 

S*(wV^W)S  =  D."^, 
J 

or 

(¥V^W)S  =   SD   "^;  (A.1+) 

t  -2 
i.e.,  S  is  the  eigenvector  matrix  of  W  A  W,  which  can  he  written  as 

W^A"^W  =  A"^  +  K,  (a. 5) 


where  the  elements  of  K  are  of  order  0(t),  i.e.,  0[|t  (y.)I/|T  (y„)|] 

m     1  m     X, 


2j 


i    ^   P.      Assume  now  A .    =   X.    ^    =    . . .    =   A,     is   an  h-i+1   fold  eigenvalue  of 

1  1+1  h 


A;    then  as   j-*^,   h-i+1   independent   eigensolutions   of   (A.i+)   with   d->X 


1 


exist.   Everyone  of  these  is  described  by  q  values  s  ,  s  ,^,  ...,  s  .   ^, 

"  p       p+1  p+q-1 

and  we  assume  that  these   solutions   are  normalized  such  that 

2  2  2 

s.      +s.,      +...+S,       =1.      Then  the   s„   with   £   ^   i,    i+1,    ...,   h  are  of 
1  1+1  h  Z 

h  P'^g"! 

order  0(t).      This  means  the  angle  between   2,   ^o   ^o    ^-"^^       I      ^o   "^o    ^^   also 

i  ^  ^        ii=p  ^  : 

h 
of  order  0(t),  while  according  to  Theorem  1  the  angle  between   J  s  w 

and  the  eigenspace  of  X.,  ^..-,>  •••»  X  of  A  is  of  order  0(q.  ).   This 

establishes  the  theorem.  | 

12 


Appendix  II 


h  = 


C C 


^8  J 


where  B  and  C  are  matrices  of  order  8, 

K. 


\  = 


a,         0.1 
k 


0.1       a,         0.1 


0.1 


0.1 
0.1      ^Na, 


,    C   =  diag(0.1+,    .  ..  ,   O.U) 


and  the   centres  of  the  Gerschgorin  discs    (a    }   are   given  "by 
(2,    3,    6,    9,    10,    11,    12,   13). 


^2  = 


16.778 
-U.889 
-1+.889 
-i;.889 
-1.556 
-0.222 


-li.889 
13.i+i+i^ 
-1.556 
-1.556 
1.778 
3.111 


-1+.889 
-1.556 

13.i+i+i+ 

-1.556 

1.778 

3.111 


-I1.889 
-1.556 
-1.556 
13.i+i+i+ 
1.778 
3.111 


-1.556 
1.778 
1.778 
1.778 

10.111 
6Mh 


-0.222 
3.111 
3.111 
3.111 
6.1+1+1+ 
8.778 


S- 


5       -h 

.1+         6 

1       -1+ 


1 

■1+         1 
6       -1+ 


-1+ 
1 


6 

-1+ 

1 


,  n  =  61+ 
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